%TQ catalog analysis
%josh crozier 2020


%% Load catalog
TQ_tab = readtable('Kilauea_2008-2018_resonant_signal_catalog_presented.csv',...
    'HeaderLines',1);

%% derivative filters
%derivative filter
Tinterp = days(1);
dfilt7 = designfilt('differentiatorfir', ... % Response type
   'FilterOrder',20, ...            % Filter order
   'PassbandFrequency',1/60, ...     % Frequency constraints
   'StopbandFrequency',1/7, ...
   'PassbandWeight',10, ...          % Design method options
   'StopbandWeight',1, ...
   'SampleRate',1/days(Tinterp));               % Sample rate (days)
% fvtool(dfilt7,'MagnitudeDisplay','zero-phase','Fs',1/days(Tinterp))
delay7 = mean(grpdelay(dfilt7));
dfilt30 = designfilt('differentiatorfir', ... % Response type
   'FilterOrder',80, ...            % Filter order
   'PassbandFrequency',1/(30*24), ...     % Frequency constraints
   'StopbandFrequency',1/40, ...
   'PassbandWeight',10, ...          % Design method options
   'StopbandWeight',1, ...
   'SampleRate',1/days(Tinterp));               % Sample rate (days)
% fvtool(dfilt30,'MagnitudeDisplay','zero-phase','Fs',1/days(Tinterp))
delay30 = mean(grpdelay(dfilt30));


%% table names for each channel
networks = "HV";
% 'using limited stations'
stationlist = ["NPB","SRM","NPT","OBL","WRM","SDH","UWE","UWB",...
    "SBL","KKO","RIMD"];
locationlist = "*";
channellist = ["HHE","HHN","HHZ"];
vertchannel_station_str = strings(1,length(stationlist));
eastchannel_station_str = strings(1,length(stationlist));
northchannel_station_str = strings(1,length(stationlist));
for statin = 1:length(stationlist)
    vertchannel_station_str(statin) = strcat(networks,'_',...
        stationlist(statin),...
        '_', channellist(3));
    eastchannel_station_str(statin) = strcat(networks,'_',...
        stationlist(statin),...
        '_', channellist(1));
    northchannel_station_str(statin) = strcat(networks,'_',...
        stationlist(statin),...
        '_', channellist(2));
end
vertchannel_station_list_amp = strcat("amp_",vertchannel_station_str);
vertchannel_station_list_phase = strcat("phase_",vertchannel_station_str);
eastchannel_station_list_amp = strcat("amp_",eastchannel_station_str);
eastchannel_station_list_phase = strcat("phase_",eastchannel_station_str);
northchannel_station_list_amp = strcat("amp_",northchannel_station_str);
northchannel_station_list_phase = strcat("phase_",northchannel_station_str);


%% mogi inversion
%
%inversion with fixed depth
TQ_tab = mogi_inv_table(TQ_tab,stationlist,...
    vertchannel_station_list_amp,vertchannel_station_list_phase,...
    eastchannel_station_list_amp, eastchannel_station_list_phase, ...
    northchannel_station_list_amp, northchannel_station_list_phase,...
    true);

%inversion with variable depth
TQ_tab = mogi_inv_table_depth(TQ_tab,stationlist,...
    vertchannel_station_list_amp,vertchannel_station_list_phase,...
    eastchannel_station_list_amp, eastchannel_station_list_phase, ...
    northchannel_station_list_amp, northchannel_station_list_phase,...
    true);

%plot example event
eventin = find(TQ_tab.start > datetime(2017,5,21,9,50,0) &...
    TQ_tab.start < datetime(2017,5,21,9,55,0) & ...
    TQ_tab.T > 34.8 & ...
    TQ_tab.T < 45);
figure()
plot_disp_map(TQ_tab(eventin,:),true,false)
%}


%% calculate mean vert and horz displacements
TQ_tab.meanvert = nanmean(TQ_tab{:,vertchannel_station_list_amp},2);
TQ_tab.meanhor = nanmean(sqrt(TQ_tab{:,eastchannel_station_list_amp}.^2 ...
+ TQ_tab{:,northchannel_station_list_amp}.^2),2);
allchannel_station_list_amp = [eastchannel_station_list_amp,...
    northchannel_station_list_amp,...
    vertchannel_station_list_amp];
TQ_tab.Amean = log10(nansum(TQ_tab{:,allchannel_station_list_amp},2)./....
    sum(isfinite(TQ_tab{:,allchannel_station_list_amp}),2));

%% calculate single station amplitude
TQ_tab.A_single_station = NaN(height(TQ_tab),1);
TQ_tab.amp_NP_E = NaN(height(TQ_tab),1);
TQ_tab.amp_NP_N = NaN(height(TQ_tab),1);
TQ_tab.amp_NP_Z = NaN(height(TQ_tab),1);
TQ_tab.phase_NP_E = NaN(height(TQ_tab),1);
TQ_tab.phase_NP_N = NaN(height(TQ_tab),1);
TQ_tab.phase_NP_Z = NaN(height(TQ_tab),1);
TQ_tab.errphase_E = NaN(height(TQ_tab),1);
TQ_tab.errphase_N = NaN(height(TQ_tab),1);
TQ_tab.errphase_Z = NaN(height(TQ_tab),1);

refdate = datetime(2008,11,1,0,0,0);
TQ_tab.A_single_station(TQ_tab.start >= refdate) = ...
    sqrt(TQ_tab.amp_HV_NPT_HHZ(TQ_tab.start >= refdate).^2);% + ...
    %TQ_tab.amp_HV_NPT_HHE(TQ_tab.start >= refdate).^2 + ...
    %TQ_tab.amp_HV_NPT_HHN(TQ_tab.start >= refdate).^2);
TQ_tab.A_single_station(TQ_tab.start < refdate) = ...
    sqrt(TQ_tab.amp_HV_SDH_HHZ(TQ_tab.start < refdate).^2);% + ...
    %TQ_tab.amp_HV_SDH_HHE(TQ_tab.start < refdate).^2 + ...
    %TQ_tab.amp_HV_SDH_HHN(TQ_tab.start < refdate).^2);

TQ_tab.amp_NP_E(TQ_tab.start >= refdate) = TQ_tab.amp_HV_NPT_HHE(TQ_tab.start >= refdate);
TQ_tab.amp_NP_N(TQ_tab.start >= refdate) = TQ_tab.amp_HV_NPT_HHN(TQ_tab.start >= refdate);
TQ_tab.amp_NP_Z(TQ_tab.start >= refdate) = TQ_tab.amp_HV_NPT_HHZ(TQ_tab.start >= refdate);
TQ_tab.amp_NP_E(TQ_tab.start < refdate) = TQ_tab.amp_HV_NPB_HHE(TQ_tab.start < refdate);
TQ_tab.amp_NP_N(TQ_tab.start < refdate) = TQ_tab.amp_HV_NPB_HHN(TQ_tab.start < refdate);
TQ_tab.amp_NP_Z(TQ_tab.start < refdate) = TQ_tab.amp_HV_NPB_HHZ(TQ_tab.start < refdate);

TQ_tab.phase_NP_E(TQ_tab.start >= refdate) = TQ_tab.phase_HV_NPT_HHE(TQ_tab.start >= refdate);
TQ_tab.phase_NP_N(TQ_tab.start >= refdate) = TQ_tab.phase_HV_NPT_HHN(TQ_tab.start >= refdate);
TQ_tab.phase_NP_Z(TQ_tab.start >= refdate) = TQ_tab.phase_HV_NPT_HHZ(TQ_tab.start >= refdate);
TQ_tab.phase_NP_E(TQ_tab.start < refdate) = TQ_tab.phase_HV_NPB_HHE(TQ_tab.start < refdate);
TQ_tab.phase_NP_N(TQ_tab.start < refdate) = TQ_tab.phase_HV_NPB_HHN(TQ_tab.start < refdate);
TQ_tab.phase_NP_Z(TQ_tab.start < refdate) = TQ_tab.phase_HV_NPB_HHZ(TQ_tab.start < refdate);

TQ_tab.errphase_E(TQ_tab.start >= refdate) = TQ_tab.errphase_HV_NPT_HHE(TQ_tab.start >= refdate);
TQ_tab.errphase_N(TQ_tab.start >= refdate) = TQ_tab.errphase_HV_NPT_HHN(TQ_tab.start >= refdate);
TQ_tab.errphase_Z(TQ_tab.start >= refdate) = TQ_tab.errphase_HV_NPT_HHZ(TQ_tab.start >= refdate);
TQ_tab.errphase_E(TQ_tab.start < refdate) = TQ_tab.errphase_HV_NPB_HHE(TQ_tab.start < refdate);
TQ_tab.errphase_N(TQ_tab.start < refdate) = TQ_tab.errphase_HV_NPB_HHN(TQ_tab.start < refdate);
TQ_tab.errphase_Z(TQ_tab.start < refdate) = TQ_tab.errphase_HV_NPB_HHZ(TQ_tab.start < refdate);

%% calculate first motions
variablenames = TQ_tab.Properties.VariableNames;
fm_z_ind = false(size(variablenames));
for ii = 1:length(variablenames)
    if length(variablenames{ii}) == 13 || length(variablenames{ii}) == 14
        if strcmp(variablenames{ii}(1:3),"fm_") && strcmp(variablenames{ii}(end-3:end),"_HHZ")
            fm_z_ind(ii) = true;
        end
    end
end
TQ_tab.fm_z_mode = mode(sign(TQ_tab{:,fm_z_ind}),2);
TQ_tab.fm_z_mode(isnan(TQ_tab.fm_z_mode)) = 0;

'!!!setting poorly determined polarity to undetermined!!!'
VLP_tab_modes.fm_z_mode(TQ_tab.fm_STALTA < 2.0) = 0;
VLP_tab_modes.fm_z_mode(TQ_tab.fm_Afrac < 0.6) = 0;


%% event density and T,Q derivatives for cr mode
% label cr mode
TQ_tab.cr = false(height(TQ_tab),1);
%T and time blocks to include
TQ_tab.cr(TQ_tab.start > datetime(2011,10,1) & ...
    TQ_tab.start < datetime(2018,7,1) & ...
    TQ_tab.T > 34.8 & TQ_tab.T < 45) = true;
TQ_tab.cr(TQ_tab.start > datetime(2011,7,7) & ...
    TQ_tab.start < datetime(2011,10,1) & ...
    TQ_tab.T > 19 & TQ_tab.T < 24) = true;
TQ_tab.cr(TQ_tab.start > datetime(2011,2,1) & ...
    TQ_tab.start < datetime(2011,3,1) & ...
    TQ_tab.T > 32 & TQ_tab.T < 37) = true;
TQ_tab.cr(TQ_tab.start > datetime(2010,9,1) & ...
    TQ_tab.start < datetime(2011,2,1) & ...
    TQ_tab.T > 28 & TQ_tab.T < 34) = true;
TQ_tab.cr(TQ_tab.start > datetime(2010,3,1) & ...
    TQ_tab.start < datetime(2010,9,1) & ...
    TQ_tab.T > 22 & TQ_tab.T < 32) = true;
TQ_tab.cr(TQ_tab.start > datetime(2009,8,23) & ...
    TQ_tab.start < datetime(2010,3,1) & ...
    TQ_tab.T > 16 & TQ_tab.T < 24) = true;
TQ_tab.cr(TQ_tab.start > datetime(2009,1,15) & ...
    TQ_tab.start < datetime(2009,8,23) & ...
    TQ_tab.T > 10 & TQ_tab.T < 20) = true;
TQ_tab.cr(TQ_tab.start > datetime(2008,7,12) & ...
    TQ_tab.start < datetime(2008,10,15) & ...
    TQ_tab.T > 20.5 & TQ_tab.T < 25.1) = true;
%blocks to exclude
TQ_tab.cr(TQ_tab.start > datetime(2012,2,2) & ...
    TQ_tab.start < datetime(2012,2,3)) = false;

%make separate table for now
keepind = TQ_tab.cr &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 2 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);
VLP_tab_cr = TQ_tab(keepind,:);
VLP_tab_cr = sortrows(VLP_tab_cr,'start');
% event densities
cr_struct = struct();
cr_struct.date = VLP_tab_cr.start(1):days(2):VLP_tab_cr.start(end);
% plotdatenum = datenum(plotdatetime);
cr_struct.density_7 = NaN(size(cr_struct.date));
cr_struct.density_30 = NaN(size(cr_struct.date));
cr_struct.T_7 = NaN(size(cr_struct.date));
cr_struct.T_30 = NaN(size(cr_struct.date));
cr_struct.Q_7 = NaN(size(cr_struct.date));
cr_struct.Q_30 = NaN(size(cr_struct.date));
for ii = 1:length(cr_struct.date)
    cr_struct.density_7(ii) = length(find(VLP_tab_cr.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(7/2)))/7;
    cr_struct.density_30(ii) = length(find(VLP_tab_cr.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(30/2)))/30;
    
    cr_struct.T_7(ii) = mean(VLP_tab_cr.T(VLP_tab_cr.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(7/2)));
    cr_struct.T_30(ii) = mean(VLP_tab_cr.T(VLP_tab_cr.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(30/2)));
    
    cr_struct.Q_7(ii) = mean(VLP_tab_cr.Q(VLP_tab_cr.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(7/2)));
    cr_struct.Q_30(ii) = mean(VLP_tab_cr.Q(VLP_tab_cr.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_cr.start <= cr_struct.date(ii) + days(30/2)));
end

cr_struct.dT_dt_7 = diff(cr_struct.T_7)./days(diff(cr_struct.date));
VLP_tab_cr.dT_dt_7 = interp1(cr_struct.date(1:end-1) + (cr_struct.date(2:end) - cr_struct.date(1:end-1))/2,...
    cr_struct.dT_dt_7, VLP_tab_cr.start);
cr_struct.dT_dt_30 = diff(cr_struct.T_30)./days(diff(cr_struct.date));
VLP_tab_cr.dT_dt_30 = interp1(cr_struct.date(1:end-1) + (cr_struct.date(2:end) - cr_struct.date(1:end-1))/2,...
    cr_struct.dT_dt_30, VLP_tab_cr.start);
cr_struct.dQ_dt_7 = diff(cr_struct.Q_7)./days(diff(cr_struct.date));
VLP_tab_cr.dQ_dt_7 = interp1(cr_struct.date(1:end-1) + (cr_struct.date(2:end) - cr_struct.date(1:end-1))/2,...
    cr_struct.dQ_dt_7, VLP_tab_cr.start);
cr_struct.dQ_dt_30 = diff(cr_struct.Q_30)./days(diff(cr_struct.date));
VLP_tab_cr.dQ_dt_30 = interp1(cr_struct.date(1:end-1) + (cr_struct.date(2:end) - cr_struct.date(1:end-1))/2,...
    cr_struct.dQ_dt_30, VLP_tab_cr.start);
VLP_tab_cr.density_7 = interp1(cr_struct.date,cr_struct.density_7,VLP_tab_cr.start);
VLP_tab_cr.density_30 = interp1(cr_struct.date,cr_struct.density_30,VLP_tab_cr.start);

%merge back with main table
VLP_tab_noncr = TQ_tab(~keepind,:);
TQ_tab = outerjoin(VLP_tab_noncr, VLP_tab_cr,'MergeKeys',true);
clear VLP_tab_noncr VLP_tab_cr

%% lake mode
% label lake mode
TQ_tab.lake = false(height(TQ_tab),1);
%T and time blocks to include
TQ_tab.lake(TQ_tab.start > datetime(2011,11,1) & ...
    TQ_tab.start < datetime(2018,7,1) & ...
    TQ_tab.T > 14 & TQ_tab.T < 20) = true;
%blocks to exclude
TQ_tab.lake(TQ_tab.start > datetime(2015,8,1) & ...
    TQ_tab.start < datetime(2018,7,1) & ...
    TQ_tab.T < 16.5) = false;

%make separate table for now
keepind = TQ_tab.lake &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 2 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);
VLP_tab_lake = TQ_tab(keepind,:);
VLP_tab_lake = sortrows(VLP_tab_lake,'start');
% event densities
lake_struct = struct();
lake_struct.date = VLP_tab_lake.start(1):days(2):VLP_tab_lake.start(end);
% plotdatenum = datenum(plotdatetime);
lake_struct.density_60 = NaN(size(lake_struct.date));
lake_struct.density_120 = NaN(size(lake_struct.date));
lake_struct.T_60 = NaN(size(lake_struct.date));
lake_struct.T_120 = NaN(size(lake_struct.date));
lake_struct.Q_60 = NaN(size(lake_struct.date));
lake_struct.Q_120 = NaN(size(lake_struct.date));
for ii = 1:length(lake_struct.date)
    lake_struct.density_60(ii) = length(find(VLP_tab_lake.start >= lake_struct.date(ii) - days(60/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(60/2)))/60;
    lake_struct.density_120(ii) = length(find(VLP_tab_lake.start >= lake_struct.date(ii) - days(120/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(120/2)))/120;
    
    lake_struct.T_60(ii) = mean(VLP_tab_lake.T(VLP_tab_lake.start >= lake_struct.date(ii) - days(60/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(60/2)));
    lake_struct.T_120(ii) = mean(VLP_tab_lake.T(VLP_tab_lake.start >= lake_struct.date(ii) - days(120/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(120/2)));
    
    lake_struct.Q_60(ii) = mean(VLP_tab_lake.Q(VLP_tab_lake.start >= lake_struct.date(ii) - days(60/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(60/2)));
    lake_struct.Q_120(ii) = mean(VLP_tab_lake.Q(VLP_tab_lake.start >= lake_struct.date(ii) - days(120/2) & ...
        VLP_tab_lake.start <= lake_struct.date(ii) + days(120/2)));
end

lake_struct.dT_dt_60 = diff(lake_struct.T_60)./days(diff(lake_struct.date));
VLP_tab_lake.dT_dt_60 = interp1(lake_struct.date(1:end-1) + (lake_struct.date(2:end) - lake_struct.date(1:end-1))/2,...
    lake_struct.dT_dt_60, VLP_tab_lake.start);
lake_struct.dT_dt_120 = diff(lake_struct.T_120)./days(diff(lake_struct.date));
VLP_tab_lake.dT_dt_120 = interp1(lake_struct.date(1:end-1) + (lake_struct.date(2:end) - lake_struct.date(1:end-1))/2,...
    lake_struct.dT_dt_120, VLP_tab_lake.start);
lake_struct.dQ_dt_60 = diff(lake_struct.Q_60)./days(diff(lake_struct.date));
VLP_tab_lake.dQ_dt_60 = interp1(lake_struct.date(1:end-1) + (lake_struct.date(2:end) - lake_struct.date(1:end-1))/2,...
    lake_struct.dQ_dt_60, VLP_tab_lake.start);
lake_struct.dQ_dt_120 = diff(lake_struct.Q_120)./days(diff(lake_struct.date));
VLP_tab_lake.dQ_dt_120 = interp1(lake_struct.date(1:end-1) + (lake_struct.date(2:end) - lake_struct.date(1:end-1))/2,...
    lake_struct.dQ_dt_120, VLP_tab_lake.start);
VLP_tab_lake.density_60 = interp1(lake_struct.date,lake_struct.density_60,VLP_tab_lake.start);
VLP_tab_lake.density_120 = interp1(lake_struct.date,lake_struct.density_120,VLP_tab_lake.start);

%merge back with main table
VLP_tab_nonlake = TQ_tab(~keepind,:);
TQ_tab = outerjoin(VLP_tab_nonlake, VLP_tab_lake,'MergeKeys',true);
clear VLP_tab_nonlake VLP_tab_lake


%% plot station timeline
%
figure('name','station timeline')
xticklabelcell = cell(size(vertchannel_station_list_amp));
for statin = 1:length(vertchannel_station_list_amp)
    plot(TQ_tab.start,statin + 0*TQ_tab{:,vertchannel_station_list_amp(statin)},...
        '-k.','MarkerSize',1), hold on
    nanplot = double(~isnan(TQ_tab{:,vertchannel_station_list_amp(statin)}));
    nanplot(logical(nanplot)) = NaN;
%     plot(TQ_tab.start,statin + nanplot,'-r.', 'MarkerSize',1), hold on
    xticklabelcell{statin} = vertchannel_station_list_amp(statin);
    xticklabelcell{statin} = xticklabelcell{statin}{1}(8:10);
end
yticks(1:length(vertchannel_station_list_amp))
yticklabels(xticklabelcell)
ylim([0.5,statin + 0.5])
%}


%% plot catalog
% keepind = TQ_tab.T > 10 & TQ_tab.T < 45 & ...
%     TQ_tab.amp_abovenoise > log10(3.0) & ...
%     TQ_tab.std_abovenoise > log10(3.0) & ...
%     TQ_tab.Q > 6 & ...
%     TQ_tab.mean_min50_errphase < 0.1 &...
%     TQ_tab.channelCount > 0 & ...
%     (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);% & ...
%     %TQ_tab.start > datetime(2011,5,1);
keepind = TQ_tab.T > 1 & TQ_tab.T < 60 & ...
    TQ_tab.amp_abovenoise > log10(2.0) & ...
    TQ_tab.std_abovenoise > log10(2.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.15 &...
    TQ_tab.channelCount > 0 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);
VLP_tab = TQ_tab(keepind,:);

figure()

numplots = 2;
posleft = 0.1;
poswidth = 0.8;
posheight = 0.8/numplots;
posheightfrac = 0.9;
posbottom = 0.1;
poscount = numplots;

Tlim = [1,50];
Xlim = [datetime(2008,1,1),datetime(2018,7,1)];
posshift = 0.04;

poscount = poscount - 1;
ax1 = subplot('Position',[posleft,posbottom+poscount*posheight + posshift,poswidth,posheight*posheightfrac]);
% fill(inflateplotx,inflateploty*100,[0.8,0.8,0.8],'EdgeColor','none')
hold on
scatter(VLP_tab.start,VLP_tab.T,4,VLP_tab.Q,'filled')
% plot(cr_struct.date,cr_struct.T_30,'k-')
grid on, grid minor
ylabel('T (s)')
cbar = colorbar;
cbar.Position = cbar.Position + [0.09,0,0,0];
cbar.Position(3) = 0.01;
colormap(ax1,cool)
caxis([6,40])
ylabel(cbar, 'Q')
% title('Conduit-Reservoir')
xlim(Xlim);%max(VLP_tab.start)])
ylim(Tlim)
xline(datetime(2008,3,19),'k')%,{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
set(gca,'Layer','top')
xticklabels([])

poscount = poscount - 1;
ax2 = subplot('Position',[posleft,posbottom+poscount*posheight + posshift,poswidth,posheight*posheightfrac]);
% fill(inflateplotx,inflateploty*100,[0.8,0.8,0.8],'EdgeColor','none')
hold on
scatter(VLP_tab.start,VLP_tab.Q,4,VLP_tab.T,'filled')
% plot(cr_struct.date,cr_struct.Q_30,'k-')
grid on, grid minor
ylabel('Q')
cbar = colorbar;
cbar.Position = cbar.Position + [0.09,0,0,0];
cbar.Position(3) = 0.01;
colormap(ax2,winter)
caxis(Tlim)
ylabel(cbar, 'T (s)')
% title('Conduit-Reservoir')
xlim(Xlim);%max(VLP_tab.start)])
ylim([6,75])
xline(datetime(2008,3,19),'k',{'Crater'},'LabelVerticalAlignment','top')
xline(datetime(2010,2,1),'k',{'SSE'},'LabelVerticalAlignment','top')
xline(datetime(2012,5,25),'k',{'SSE'},'LabelVerticalAlignment','top')
xline(datetime(2011,3,5),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2011,8,3),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2011,9,21),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2014,6,27),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2016,5,24),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2018,4,30),'k',{'ERZ'},'LabelVerticalAlignment','top')
xline(datetime(2012,10,26),'k',{'Int'},'LabelVerticalAlignment','top')
xline(datetime(2014,5,10),'k',{'Int'},'LabelVerticalAlignment','top')
xline(datetime(2015,5,9),'k',{'Int'},'LabelVerticalAlignment','top')
set(gca,'Layer','top')
% xticklabels([])


%% plot conduit-reservoir and lake modes
keepind = TQ_tab.T > 10 & TQ_tab.T < 45 &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 0 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2) & ...
    TQ_tab.start > datetime(2011,10,1);
VLP_tab_modes = TQ_tab(keepind,:);

% VLP_tab_modes.fm_z_mode(VLP_tab_modes.T < 34.8) = 0;
VLP_tab_modes.fm_z_mode(~VLP_tab_modes.cr) = 0;
VLP_tab_modes.fm_z_mode(TQ_tab.fm_STALTA < 2.5) = 0;
VLP_tab_modes.fm_z_mode(TQ_tab.fm_Afrac < 0.6) = 0;

VLP_tab_modes.isolated = true(height(VLP_tab_modes),1);
VLP_tab_modes.label = strings(height(VLP_tab_modes),1);
VLP_tab_modes.trigger = VLP_tab_modes.fm_z_mode;
cr_amp = [];
cr_trigger = [];
cr_Q = [];
lake1_amp = [];
lake2_amp = [];
lake1_T = [];
lake1_Q = [];
lake1_crQ = [];
lake2_T = [];
timegap = seconds(120);
startin = 1;
VLP_tab_modes = sortrows(VLP_tab_modes,'start');
event_groups = {};
for ii = 1:height(VLP_tab_modes)
    if abs(VLP_tab_modes.start(ii) - VLP_tab_modes.start(startin)) > timegap
        if ii - startin > 1
            ivec_temp = startin:ii-1;
            VLP_tab_modes.isolated(ivec_temp) = false;
            event_groups{end+1} = ivec_temp;
            group_T = VLP_tab_modes.T(event_groups{end});
            %cr = find(group_T > 34.8,1);
            cr = find(VLP_tab_modes.cr(event_groups{end}));
            lake_in = find(group_T > 12 & group_T < 22);
            if length(lake_in) == 1
                lake1 = lake_in;
                lake2 = [];
            elseif length(lake_in) > 1
                [~,lake1_in] = max(group_T(lake_in));
                lake1 = lake_in(lake1_in);
                lake_in_trim = lake_in([1:lake1_in-1,lake1_in+1:end]);
                [~,lake2_in] = max(group_T(lake_in_trim));
                lake2 = lake_in_trim(lake2_in);
            else
                lake1 = [];
                lake2 = [];
            end
            if ~isempty(cr) && ~isempty(lake1)
                cr_amp(end+1) = VLP_tab_modes.A_single_station(event_groups{end}(cr));
                cr_Q(end+1) = VLP_tab_modes.Q(event_groups{end}(cr));
                VLP_tab_modes.label(event_groups{end}(cr)) = "cr";
                cr_trigger(end+1) = VLP_tab_modes.fm_z_mode(event_groups{end}(cr));
                VLP_tab_modes.trigger(event_groups{end}(cr)) = VLP_tab_modes.fm_z_mode(event_groups{end}(cr));
                lake1_amp(end+1) = VLP_tab_modes.A_single_station(event_groups{end}(lake1));
                lake1_T(end+1) = VLP_tab_modes.T(event_groups{end}(lake1));
                lake1_Q(end+1) = VLP_tab_modes.Q(event_groups{end}(lake1));
                lake1_crQ(end+1) = VLP_tab_modes.Q(event_groups{end}(cr));
                VLP_tab_modes.label(event_groups{end}(lake1)) = "lake1";
                VLP_tab_modes.trigger(event_groups{end}(lake1)) = VLP_tab_modes.fm_z_mode(event_groups{end}(cr));
                if ~isempty(lake2)
                    lake2_amp(end+1) = VLP_tab_modes.A_single_station(event_groups{end}(lake2));
                    lake2_T(end+1) = VLP_tab_modes.T(event_groups{end}(lake2));
                    VLP_tab_modes.label(event_groups{end}(lake2)) = "lake2";
                    VLP_tab_modes.trigger(event_groups{end}(lake2)) = VLP_tab_modes.fm_z_mode(event_groups{end}(cr));
                else
                    lake2_amp(end+1) = NaN;
                    lake2_T(end+1) = NaN;
                end
            end
        end
        startin = ii;
    end
end

cr_struct.density_7 = NaN(size(cr_struct.date));
cr_struct.density_30 = NaN(size(cr_struct.date));
cr_struct.densitynormal_7 = NaN(size(cr_struct.date));
cr_struct.densitynormal_30 = NaN(size(cr_struct.date));
cr_struct.densityreverse_7 = NaN(size(cr_struct.date));
cr_struct.densityreverse_30 = NaN(size(cr_struct.date));
for ii = 1:length(cr_struct.date)
    cr_struct.density_7(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(7/2) & ...
        VLP_tab_modes.cr))/7;
    cr_struct.density_30(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(30/2) & ...
        VLP_tab_modes.cr))/30;
    
    cr_struct.densitynormal_7(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(7/2) & ...
        VLP_tab_modes.fm_z_mode == 1))/7;
    cr_struct.densitynormal_30(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(30/2) & ...
        VLP_tab_modes.fm_z_mode == 1))/30;
    
    cr_struct.densityreverse_7(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(7/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(7/2) & ...
        VLP_tab_modes.fm_z_mode == -1))/7;
    cr_struct.densityreverse_30(ii) = length(find(VLP_tab_modes.start >= cr_struct.date(ii) - days(30/2) & ...
        VLP_tab_modes.start <= cr_struct.date(ii) + days(30/2) & ...
        VLP_tab_modes.fm_z_mode == -1))/30;
end

figure()
subplot(3,2,[1,2])
p1 = plot(VLP_tab_modes.start(VLP_tab_modes.cr & VLP_tab_modes.trigger == 0),...
    VLP_tab_modes.T(VLP_tab_modes.cr & VLP_tab_modes.trigger == 0),'k.');
hold on
p2 = plot(VLP_tab_modes.start(VLP_tab_modes.cr & VLP_tab_modes.trigger == 1),...
    VLP_tab_modes.T(VLP_tab_modes.cr & VLP_tab_modes.trigger == 1),'g.');
p3 = plot(VLP_tab_modes.start(VLP_tab_modes.cr & VLP_tab_modes.trigger == -1),...
    VLP_tab_modes.T(VLP_tab_modes.cr & VLP_tab_modes.trigger == -1),'.','MarkerEdgeColor',[0.8500, 0.3250, 0.0980]);
p5 = plot(VLP_tab_modes.start(~VLP_tab_modes.cr),...
    VLP_tab_modes.T(~VLP_tab_modes.cr),'ko','MarkerSize',3);
p4 = plot(VLP_tab_modes.start(strcmp(VLP_tab_modes.label,"lake1")),...
    VLP_tab_modes.T(strcmp(VLP_tab_modes.label,"lake1")),'go','MarkerSize',3);
ylabel('T (s)')
legend([p2,p3,p1,p4,p5],{'Normal CR Mode','Reverse CR Mode','Undetermined CR Mode',...
    'Concurrent Lake Mode','Other/Isolated Mode'},'AutoUpdate','off')
grid on
xline(datetime(2008,3,19),'k',{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k',{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k',{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k',{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k',{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k',{'Int'},'LabelVerticalAlignment','bottom')
xlim([datetime(2011,9,1),datetime(2018,7,1)])
set(gca,'Layer','top')
ylim([10,43])

subplot(3,2,[3,4])
p2 = plot(cr_struct.date,cr_struct.density_30,'-k','LineWidth',2);
hold on
p2N = plot(cr_struct.date,cr_struct.densitynormal_30,'-g','LineWidth',1);
p2R = plot(cr_struct.date,cr_struct.densityreverse_30,'-','LineWidth',1,'Color',[0.8500, 0.3250, 0.0980]);
ylabel('Events/day')
legend([p2,p2N,p2R],["Total Conduit-Reservoir events/day",...
    "Normal events/day","Reverse events/day"],'Autoupdate','off')
xline(datetime(2008,3,19),'k')%,{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xlim([datetime(2011,9,1),datetime(2018,7,1)])
ylim([0, max(cr_struct.density_30)])
set(gca,'Layer','top')
grid on

ax2 = subplot(3,2,5);
cr_trigger(isnan(cr_trigger)) = 0;
p2 = scatter(cr_amp(cr_trigger == 1),lake1_amp(cr_trigger == 1),9,lake1_Q(cr_trigger == 1),'filled');%lake1_crQ(cr_trigger == 1),'filled');
hold on
xlabel('Conduit-Reservoir NPT Z Amp (m/s)')
ylabel('Lake Sloshing NPT Z Amp (m/s)')
set(gca,'XScale','log')
set(gca,'YScale','log')
cbar = colorbar;
ylabel(cbar,'Lake Sloshing Q')
colormap(ax2,cool)
grid on
xlim([min(cr_amp),max(cr_amp)])
ylim([min(lake1_amp),max(lake1_amp)])

ax3 = subplot(3,2,6);
cr_trigger(isnan(cr_trigger)) = 0;
p2 = scatter(cr_Q(cr_trigger == 1),lake1_Q(cr_trigger == 1),9,log10(lake1_amp(cr_trigger == 1)),'filled');%lake1_crQ(cr_trigger == 1),'filled');
hold on
xlabel('Conduit-Reservoir Q')
ylabel('Lake Sloshing Q')
cbar = colorbar;
ylabel(cbar,'Lake Sloshing NPT Z Amp')
colormap(ax3,summer)
grid on
% xlim([min(cr_Q),max(cr_Q)])
ylim([min(lake1_Q),max(lake1_Q)])



%% correlation
keepind = TQ_tab.cr &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 0 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2) & ...
    TQ_tab.start > datetime(2011,10,1);% &...TQ_tab.fm_z_mode == -1;
% keepind = TQ_tab.cr &...
%     TQ_tab.amp_abovenoise > log10(3.0) & ...
%     TQ_tab.std_abovenoise > log10(3.0) & ...
%     TQ_tab.Q > 6 & ...
%     TQ_tab.mean_min50_errphase < 0.1 &...
%     TQ_tab.channelCount > 2 & ...
%     (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);
VLP_tab_cr = TQ_tab(keepind,:);
VLP_tab_cr = sortrows(VLP_tab_cr,'start');

VLP_tab_cr.verthor = log10(VLP_tab_cr.meanvert./VLP_tab_cr.meanhor);
corrtable = VLP_tab_cr(:,{'T','Q','dT_dt_7','dQ_dt_7'...
    'A_single_station',...
    'mogidepth_depth','mogi_mean_horangle_misfit'});
corrtable.A_single_station = log10(corrtable.A_single_station);
corrtable.mogidepth_depth = corrtable.mogidepth_depth/10^3;
figure()
corrplot_mod(corrtable,'testR','on','alpha',0.05,'type','Pearson',...
    'varNames',{'T \newline(s)','Q \newline ', 'dT/dt \newline(s/day)', 'dQ/dt \newline(/day)'...
    'NPT Z Amp \newline(log_{10} m/s)',...
    'Mogi Depth \newline(km)','Radial Misfit\newline(deg)'});

%% moving correlation
keepind = TQ_tab.cr &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 0 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);
VLP_tab_cr = TQ_tab(keepind,:);
VLP_tab_cr = sortrows(VLP_tab_cr,'start');

movcorr.Tinterp = days(2);
tspan = 90;
movcorr.date_interp = datetime(2008,1,1):movcorr.Tinterp:datetime(2018,6,1);
movcorr.T_Q_30 = NaN(size(movcorr.date_interp));
for ii = 1:length(movcorr.date_interp)
    Tdata = VLP_tab_cr.T(VLP_tab_cr.start < movcorr.date_interp(ii) + days(tspan/2) & ...
        VLP_tab_cr.start > movcorr.date_interp(ii) - days(tspan/2)); 
    Qdata = VLP_tab_cr.Q(VLP_tab_cr.start < movcorr.date_interp(ii) + days(tspan/2) & ...
        VLP_tab_cr.start > movcorr.date_interp(ii) - days(tspan/2)); 
    
    if length(Tdata) > 2
        [movcorr.T_Q_30(ii),movcorr.T_Q_30_p(ii)] = corr(Tdata,Qdata);
    end
end

figure()
scatter(movcorr.date_interp,movcorr.T_Q_30,(1 - movcorr.T_Q_30_p)*6,movcorr.T_Q_30,'filled')
hold on
colormap(cmocean('balance'))
caxis([-1,1])
% cbar = colorbar;
% ylabel(cbar,'P Value')
grid on
set(gca,'layer','top')
ylim([-1,1])
xlim([datetime(2008,1,1),datetime(2018,7,1)])
xline(datetime(2008,3,19),'k')%,{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
title('T vs Q')
ylabel('Correlation \newline Coefficient')
yline(0,'k')


%% normal vs reverse statistics
keepind = TQ_tab.cr &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 0 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2) & ...
    TQ_tab.start > datetime(2011,10,1) &...
    TQ_tab.fm_z_mode ~= 0;
VLP_tab_cr = TQ_tab(keepind,:);
VLP_tab_cr = sortrows(VLP_tab_cr,'start');

figure()
%edges = linspace(min(VLP_tab_cr.T), max(VLP_tab_cr.T),15);
edges = [sort(unique(VLP_tab_cr.T)-.001); max(VLP_tab_cr.T) + .001];
subplot(2,8,1)
histogram(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == 1), edges,...
    'EdgeColor','none', 'FaceColor','k')
p1 = xline(nanmean(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == 1)),'b','LineWidth',1.5);
p2 = xline(nanmedian(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == 1)),'r:','LineWidth',1.5);
% legend([p1,p2],{'Mean','Median'})
grid on
ylabel('   Normal \newline Conduit-Reservoir Event \newline Count (2012-2018)')
subplot(2,8,9)
histogram(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == -1), edges,...
    'EdgeColor','none', 'FaceColor','k')
xline(nanmean(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == -1)),'b','LineWidth',1.5)
xline(nanmedian(VLP_tab_cr.T(VLP_tab_cr.fm_z_mode == -1)),'r:','LineWidth',1.5)
grid on
ylabel('   Reverse \newline Conduit-Reservoir Event \newline Count (2012-2018)')
xlabel('T\newline(s)')

edges = linspace(min(VLP_tab_cr.Q), 70,20);
subplot(2,8,2)
histogram(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == 1), edges,...
    'EdgeColor','none', 'FaceColor','k')
p1 = xline(nanmean(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == 1)),'b','LineWidth',1.5);
p2 = xline(nanmedian(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == 1)),'r:','LineWidth',1.5);
legend([p1,p2],{'Mean','Median'})
grid on
% ylabel('Normal Conduit-Reservoir Event \newline Count (2012-2018)')
subplot(2,8,10)
histogram(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == -1), edges,...
    'EdgeColor','none', 'FaceColor','k')
xline(nanmean(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == -1)),'b','LineWidth',1.5)
xline(nanmedian(VLP_tab_cr.Q(VLP_tab_cr.fm_z_mode == -1)),'r:','LineWidth',1.5)
grid on
% ylabel('Reverse Conduit-Reservoir Event \newline Count (2012-2018)')
xlabel('Q')

edges = linspace(log10(min(VLP_tab_cr.A_single_station)), log10(max(VLP_tab_cr.A_single_station)),20);
subplot(2,8,3)
histogram(log10(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == 1)), edges,...
    'EdgeColor','none', 'FaceColor','k')
xline(log10(nanmean(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == 1))),'b','LineWidth',1.5)
xline(nanmedian(log10(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == 1))),'r:','LineWidth',1.5)
grid on
%ylabel('Normal Count')
subplot(2,8,11)
histogram(log10(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == -1)), edges,...
    'EdgeColor','none', 'FaceColor','k')
xline(log10(nanmean(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == -1))),'b','LineWidth',1.5)
xline(nanmedian(log10(VLP_tab_cr.A_single_station(VLP_tab_cr.fm_z_mode == -1))),'r:','LineWidth',1.5)
grid on
%ylabel('Reverse Count')
xlabel('log_{10} NPT Z Amplitude\newline (m/s)')



%% plot displacement patterns

keepind = TQ_tab.cr &...
    TQ_tab.amp_abovenoise > log10(3.0) & ...
    TQ_tab.std_abovenoise > log10(3.0) & ...
    TQ_tab.Q > 6 & ...
    TQ_tab.mean_min50_errphase < 0.1 &...
    TQ_tab.channelCount > 2 & ...
    (TQ_tab.Acwt - TQ_tab.maxamp) > log10(2);% & ...TQ_tab.start > datetime(2012,6,1);

'!!!excluding bad timeblock from movavg plot!!!'
keepind(TQ_tab.start > datetime(2011,10,1) & ...
    TQ_tab.start < datetime(2012,7,1)) = false;


keepind(TQ_tab.start > datetime(2010,10,1) & ...
    TQ_tab.channelCount < 18) = false;
keepind(TQ_tab.channelCount > 5 & ...
    TQ_tab.channelCount < 18) = false;


VLP_tab_cr = TQ_tab(keepind,:);
VLP_tab_cr = sortrows(VLP_tab_cr,'start');

ind1 = VLP_tab_cr.channelCount < 4;
ind2 = VLP_tab_cr.channelCount >= 4 & VLP_tab_cr.channelCount < 18;
ind3 = VLP_tab_cr.channelCount >= 18;

plotdatetime = VLP_tab_cr.start(1):days(2):VLP_tab_cr.start(end);

plotvar = VLP_tab_cr.mogidepth_depth;
colorvar = log10(VLP_tab_cr.mogidepth_misfit*100);
density_120_ind1 = NaN(size(plotdatetime));
density_120_ind3 = NaN(size(plotdatetime));
for ii = 1:length(plotdatetime)
    
    dataR = plotvar(VLP_tab_cr.start >= plotdatetime(ii) - days(120/2) & ...
        VLP_tab_cr.start <= plotdatetime(ii) & ind1);
    dataL = plotvar(VLP_tab_cr.start <= plotdatetime(ii) + days(120/2) & ...
         VLP_tab_cr.start >= plotdatetime(ii) & ind1);
    if ~(isempty(dataR) || isempty(dataL))
        density_120_ind1(ii) = mean([dataL;dataR]);
    end
    dataR = plotvar(VLP_tab_cr.start >= plotdatetime(ii) - days(120/2) & ...
        VLP_tab_cr.start <= plotdatetime(ii) & ind3);
    dataL = plotvar(VLP_tab_cr.start <= plotdatetime(ii) + days(120/2) & ...
         VLP_tab_cr.start >= plotdatetime(ii) & ind3);
    if ~(isempty(dataR) || isempty(dataL))
        density_120_ind3(ii) = mean([dataL;dataR]);
    end
end

figure()
ax2 = subplot(2,1,2);
scatter(VLP_tab_cr.start(ind1), plotvar(ind1),6,colorvar(ind1),'x')
hold on
set(gca,'layer','top')
scatter(VLP_tab_cr.start(ind3), plotvar(ind3),6,colorvar(ind3),'filled')
hold on
plot(plotdatetime,density_120_ind1,'r','LineWidth',1)
plot(plotdatetime,density_120_ind3,'k','LineWidth',1)
ylabel("Mogi Depth (m)")
grid on
grid minor
xline(datetime(2008,3,19),'k',{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k',{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k',{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k',{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k',{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k',{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k',{'Int'},'LabelVerticalAlignment','bottom')
cbar = colorbar;
ylabel(cbar,'log_{10} Mogi Misfit %')
colormap(ax2,cmocean('thermal'))
caxis([log10(10),log10(120)])
ylim([500,2500])
xlim([datetime(2008,1,1),datetime(2018,7,1)])
title('Conduit-Reservoir Mode Mogi Inverions')

ax1 = subplot(2,1,1);
plotvar = log10(VLP_tab_cr.meanvert./VLP_tab_cr.meanhor);
colorvar = log10(VLP_tab_cr.mogi_mean_horangle_misfit);
density_120_ind1 = NaN(size(plotdatetime));
density_120_ind3 = NaN(size(plotdatetime));
for ii = 1:length(plotdatetime)
    
    dataR = plotvar(VLP_tab_cr.start >= plotdatetime(ii) - days(120/2) & ...
        VLP_tab_cr.start <= plotdatetime(ii) & ind1);
    dataL = plotvar(VLP_tab_cr.start <= plotdatetime(ii) + days(120/2) & ...
         VLP_tab_cr.start >= plotdatetime(ii) & ind1);
    if ~(isempty(dataR) || isempty(dataL))
        density_120_ind1(ii) = mean([dataL;dataR]);
    end
    dataR = plotvar(VLP_tab_cr.start >= plotdatetime(ii) - days(120/2) & ...
        VLP_tab_cr.start <= plotdatetime(ii) & ind3);
    dataL = plotvar(VLP_tab_cr.start <= plotdatetime(ii) + days(120/2) & ...
         VLP_tab_cr.start >= plotdatetime(ii) & ind3);
    if ~(isempty(dataR) || isempty(dataL))
        density_120_ind3(ii) = mean([dataL;dataR]);
    end
end


scatter(VLP_tab_cr.start(ind1), plotvar(ind1),6,colorvar(ind1),'x')
hold on
set(gca,'layer','top')
scatter(VLP_tab_cr.start(ind3), plotvar(ind3),6,colorvar(ind3),'filled')
hold on
plot(plotdatetime,density_120_ind1,'r','LineWidth',1)
plot(plotdatetime,density_120_ind3,'k','LineWidth',1)
ylabel("log_{10} Vertical/Horizontal")
grid on
grid minor
xline(datetime(2008,3,19),'k')%,{'Crater'},'LabelVerticalAlignment','bottom')
xline(datetime(2010,2,1),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,5,25),'k')%,{'SSE'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,3,5),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,8,3),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2011,9,21),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,6,27),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2016,5,24),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2018,4,30),'k')%,{'ERZ'},'LabelVerticalAlignment','bottom')
xline(datetime(2012,10,26),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2014,5,10),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
xline(datetime(2015,5,9),'k')%,{'Int'},'LabelVerticalAlignment','bottom')
cbar = colorbar;
ylabel(cbar,'log_{10} Radial Misfit (degrees)')
colormap(ax1,cmocean('thermal'))
caxis(log10([5,60]))
ylim([-1,max(plotvar)])
xlim([datetime(2008,1,1),datetime(2018,7,1)])
title('Conduit-Reservoir Mode Displacement Patterns')
xticklabels([])


